Collective behavior of interacting 
self-propelled particles 



Andras Czirok and Tamas Vicsek 

Dept. of Biological Physics, Eotvos University, Budapest, Hungary 



Abstract 

We discuss biologically inspired, inherently non-equilibrium self-propelled particle 
models, in which the particles interact with their neighbours by choosing at each 
time step the local average direction of motion. We summarize some of the results 
of large scale simulations and theoretical approaches to the problem. 



1 Introduction 

The collective motion of various organisms like the flocking of birds, swimming 
of schools of fish, motion of herds of quadrupeds (see, e.g. [1,2] and references 
therein) migrating bacteria [3], molds [4], ants [5] or pedestrians [6] is a fas- 
cinating phenomenon of nature. We address the question whether there are 
some global, perhaps universal features of this type of behavior when many 
organisms are involved and parameters like the level of perturbations or the 
mean distance between the individuals are changed. 

These studies are also motivated by recent developments in areas related to 
statistical physics. Concepts originated from the physics of phase transitions 
in equilibrium systems [7] such as scale invariance and renormalization have 
also been shown to be useful in the understanding of various non-equilibrium 
systems, typical in our natural and social environment. Motion and related 
transport phenomena represent further characteristic aspects of many non- 
equilibrium processes and they are essential features of most living systems. 

To study the collective motion of large groups of organisms, the concept of 
self propelled particle (SPP) models was introduced in [8]. As the motion of 
flocking organisms is usually controlled by interactions with their neighbors 
[1], the SPP models consist of locally interacting particles with an intrinsic 
driving force, hence with a finite steady velocity. Because of their simplic- 
ity, such models represent a statistical approach complementing other studies 
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which take into account more details of the actual behavior [6,9,10], but treat 
only a moderate number of organisms and concentrate less on the large scale 
behavior. 

In spite of the analogies with ferromagnetic models, the general behavior of 
SPP systems can be quite different from those observed in equilibrium mag- 
nets. In particular, equilibrium ferromagnets possessing continuous rotational 
symmetry do not have ordered phase at finite temperatures in two dimensions 
[11]. However, in 2d SPP models an ordered phase can exist at finite noise 
levels (temperatures) as it was first demonstrated by simulations [8,12] and 
explained by a theory of flocking developed by Toner and Tu [13,14]. Further 
studies revealed that modeling collective motion leads to interesting specific 
results in all of the relevant dimensions (from 1 to 3) [15,16]. 



2 Discrete models 

The 2d system. The simplest model, introduced in [8], consists of particles 
moving on a plane with periodic boundary condition. The particles are char- 
acterized by their (off-lattice) location X{ and velocity V{ pointing in direction 
The self-propelled nature of the particles is manifested by keeping the mag- 
nitude of the velocity fixed to vq. Particles interact through the following local 
rule: at each time step a given particle assumes the average direction of motion 
of the particles in its local neighborhood S(i) (e.g., in a circle of some given 
radius centered at the position of the ith particle) with some uncertainty, as 
described by 



where the noise £ is a random variable with a uniform distribution in the 
interval [—77/2,77/2]. The locations of the particles are updated as 



with \vi\ — v = const. 

The model defined by Eqs. (1) and (2) is a transport related, non-equilibrium 
analog of ferromagnetic models [17]. The analogy is as follows: the Hamilto- 
nian tending to align the spins in the same direction in the case of equilibrium 
ferromagnets is replaced by the rule of aligning the direction of motion of 
particles, and the amplitude of the random perturbations can be considered 
proportional to the temperature for 77 <C 1. From a hydrodynamical point of 
view, in SPP systems the momentum of the particles is not conserved and - 




Xi(t + At) 



Xi(t)+Vi(t)At 



(2) 



2 



as was pointed out in [18] - the Galilean invariance is broken. Thus, the flow 
fields emerging in these models can considerably differ from the usual behavior 
of fluids. 



Critical exponents. The model defined through Eqs. (1) and (2) with cir- 
cular interaction range was studied by performing large-scale simulations. Due 
to the simplicity of the model, only two control parameters should be distin- 
guished: the (average) density of particles g = N/ L 2 (where N is the number 
of particles and L is the system size in units of the interaction range) and 
the amplitude of the noise r\. Depending on the value of these parameters the 
model can exhibit various type of behaviors as Fig. 1 demonstrates. 

For the statistical characterization of the system a well-suited order parameter 
is the magnitude of the average momentum of the particles: <fi = J2j Vj . 
This measure of the net flow is non-zero in the ordered phase, and vanishes (for 
an infinite system) in the disordered phase. Since the simulations were started 
from a random, disordered configuration, <p(t = 0) ~ 0. After some relaxation 
time a steady state emerges indicated by the convergence of the cumulative 
average (1/r) J T 4>(i)dt. The stationary values of are plotted in Fig. 2 vs r\ for 
various parameters. For weak noise the model displays an ordered motion, i.e. 
</> ~ 1, which disappears in a continuous manner by increasing r\. As L — > oo, 
the numerical results show the presence of a phase transition described by 



where r) c (g) is the critical noise amplitude that separates the ordered and 
disordered phases and (3 2 d = 0.42 ± 0.03, was found [12] to be different from 
the the mean-field value 1/2. This large-scale behavior does not depend on 
the specific choice of the microscopic details as Fig. 2 demonstrates. 

As a further analogy with equilibrium phase transitions, the fluctuations of 
the order parameter also increase on approaching the critical line [12]. The 
tails of the curves are symmetric, and decay as power-laws with an exponent 
7 close to 2, which value is, again, different from the mean-field result. 

Next we discuss the role of density. The long-range ordered phase is present 
for any g, but for a fixed value of r], vanishes with decreasing g. The critical 
line rj c (g) in the 77 — g parameter space was found to follow 



with k = 0.45 ± 0.05 for p < 1. This critical line is qualitatively different 
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vdg) ~ g K , 
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Fig. 1. Typical configurations of the positions of SPPs displayed for various values 
of density and noise. For small enough noise coherently moving clusters can be iden- 
tified in the system, which gradually disappear as the noise amplitude is increased. 
The noise amplitudes applied in the simulations were the following: n = 0.2 (a,d); 
rj = 0.4 (b,e); n = 0.8 (c,f); n = 0.6 (g); n = 1.2 (h) and n = 1.8 (i). The corre- 
sponding order parameter values are: <p = 0.75 (a,d); (ft = 0.19 (b); (ft = 0.025 (c); 
(ft = 0.3 (e,h); (ft = 0.05 (f,i) and (ft = 0.85 (g). The scales are given in units of the 
interaction length. 



from that of the diluted ferromagnets, since here the critical density at r] — > 
(corresponding to the percolation threshold for diluted ferromagnets, see, e.g., 
[17]) is vanishing. 

These findings indicate that SPP systems can be quite well characterized us- 
ing the framework of classical critical phenomena, but also show surprising 
features when compared to the analogous equilibrium systems. The velocity 
Vq provides a control parameter which switches between the SPP behavior 
(i>o > 0) and an XY ferromagnet (t>o = 0). Indeed, for v q = Kosterlitz- 
Thouless vortices [19] can be observed in the system, which are unstable 
(Fig. 3.) for any nonzero t> investigated in [12]. 
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Fig. 2. Left: Schematic illustration of the various realizations of the SPP model. 
Particles move off-lattice on a plane and interact with other particles located in the 
surrounding, which can be either a circle (Si) or 9 neighboring cells in an underlying 
lattice (Si). We plot these interaction areas for particle A with a solid and dashed 
line, respectively. Right: The average momentum of the system in the steady state vs 
the rescaled noise amplitude n/n c for various interaction ranges, densities and system 
sizes. Open symbols refer to interaction type S± and q = 1 while filled symbols refer 
to interaction type S2 and q = 2. The order present at small rj disappears in a 
continuous manner reminiscent of second order phase transitions for N — > 00. 

Correlation functions. Beside the calculation of the order parameter, it 
is insightful to characterize the configurations with correlation functions, such 
as the velocity-velocity correlation function 

C (r) = (v{r + 7,t)v{7,t))?,j (5) 



where v(f, t) is the coarse-grained velocity field and the average is taken over 
all possible values of f and t. The system is ordered on the macroscopic 
scale if = lim^^oo C(r) > 0. The decay of the fluctuations is given by 
the connected piece of the correlation function Cc(f), defined as Ccif) = 

One of the major result of the analysis of Toner and Tu [14] was that the 
attenuation of fluctuations is anisotropic in a SPP flock. In particular, they 
predicted that within a certain length scale \r\ < £, C(r) decays as 

C(f) ~ r^ 2/5 (6) 



with r± and rn being the orthogonal and parallel projection of f relative to 
the average direction of motion {v{f,t)) ?l , respectively. 

We calculated the equal-time velocity- velocity correlation function in (r±, rn) 
base for various time moments and averaged them to obtain C(r). To demon- 
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Fig. 3. Snapshots of the time development of a system with N = 4000, L = 40 
and v = 0.01 at 50 (a), 100 (b), 400 (c) and 3000 (d) steps. First the behavior is 
reminiscent of the equilibrium XY model, where the long range order is missing 
since vortices are present in the system. Here the vortices are unstable, and finally 
a self-organized long-range order develops. (After [12].) 

strate the predicted anisotropy of the correlation functions, in Fig. 4 curves 
are plotted which represent averages of points for which either the r± > r\\ or 
the r± < r| | relation holds: 

Ci.tr) = (C(f)) |r1=rjr±>r|| , q,(r) = (C(f)) |r1=rjr±<r|| . (7) 

For r] < T) c the anisotropy of the correlation functions is clearly seen, while the 
behavior of C± is consistent with the predicted (6). In contrast, for rj > rj c the 
curves vanish for large r and the system is isotropic as expected. 

Role of boundary conditions. Simulations [20-22] with reflective bound- 
aries and a short range repulsive force (constraining the maximal local density 
in the system) pointed to the importance of the boundary conditions in SPP 
models. 
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Fig. 4. Equal time velocity-velocity correlation functions characterizing the decay 
of correlations parallel and perpendicular to the average direction of motion of 
the hock. C\\ and C± were calculated for systems with circular interaction range, 
N = 4000, L = 480 and n = 0.02 (a), rj = 0.32 (b). At these parameter values the 
critical noise amplitude is at r\ c k, 0.2. The solid line in (a) represents a power-law 
ht on Cc with the predicted exponent —2/5. 




Fig. 5. A possible stationary state of the model with reflective boundary conditions. 
(After [22]) 

In such simulations rotation of the particles develop (see Fig. 2) in the high 
density, low noise regime. The direction of the rotation is selected by spon- 
taneous symmetry breaking, thus both clockwise and anti-clockwise spinning 
"vortices" can emerge. This rotating state should be distinguished from the 
Kosterlitz-Thouless vortices [19], since in our case a single vortex develops 
irrespective to the system size. In Ref. [22] we give examples of such vortices 
developing in nature. 



One and three dimensional SPP systems. Since in Id the particles 
cannot get around each other, some of the important features of the dynamics 
present in higher dimensions are lost. On the other hand, motion in Id implies 
new interesting aspects (groups of the particles have to be able to change their 
direction for the opposite in an organized manner) and the algorithms used 
for higher dimensions should be modified to take into account the specific 
crowding effects typical for Id (the particles can slow down before changing 
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direction and dense regions may be built up of momentarily oppositely moving 
particles). 

In a way the system studied below can be considered as a model of people, 
moving in a narrow channel. Imagine that a fire alarm goes on, the tunnel is 
dark, smoky, everyone is extremely excited. People are both trying to follow 
the others (to escape together) and behave in an erratic manner (due to smoke 
and excitement). 

Thus, in [15] N off-lattice particles along a line of length L have been consid- 
ered. The particles are characterized by their coordinate Xj and dimensionless 
velocity Ui updated as 



Xi(t + At) = Xi (t) + v Ui(t)At, (8) 
Ul (t + At) =G((u(t)) s( i)) + (9) 

The local average velocity (u)su) for the ith particle is calculated over the 
particles located in the interval [xj — A, Xj+A], where A = 1/2. The function G 
incorporates both the "propulsion" and "friction" forces which set the velocity 
to a prescribed value vq on the average: G(u) > u for u < 1 and G(u) < u for 
u > 1. In the numerical simulations [15] one of the simplest choices for G was 
implemented as 

v 1 I (u - l)/2 for u < 0, v 1 



and random initial and periodic boundary conditions were applied. 

Again, the emergence of the ordered phase was observed through a second 
order phase transition [15] with (3\d = 0.60 ±0.05, which is different from both 
the the mean-field value 1/2 and fad ~ 0.4 found in 2d. The critical line on 
the p — 7] phase diagram also follows (4) with re ld 1/4. 

In two dimensions an effective long range interaction can build up because 
the migrating particles have a considerably higher chance to get close to each 
other and interact than in three dimensions (where, as is well known, random 
trajectories do not overlap). The less interaction acts against ordering. On the 
other hand, in three dimensions even regular ferromagnets order. Thus, it is 
interesting to see how these two competing features change the behavior of 3d 
SPP systems. The convenient generalization of (1) for the 3d case can be the 
following [16]: 

Vi(t + At) = v N( N((v(t)) s{i) ) + £), (11) 
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where N(tt) = u/\u\ and the noise £ is uniformly distributed in a sphere of 
radius rj. 

Generally, the behavior of the system was found [16] to be similar to that of 
described above. The long-range ordered phase was present for any g, but for 
a fixed value of rj, <fi vanished with decreasing g. To compare this behavior 
to the corresponding diluted ferromagnet, <f)(r], g) was determined for t> — 0, 
when the model reduces to an equilibrium system of randomly distributed 
"spins" with a ferromagnetic-like interaction. Again, a major difference was 
found between the SPP and the equilibrium models: in the static case the 
system does not order for densities below a critical value close to 1 which 
corresponds to the percolation threshold of randomly distributed spheres in 
3d. 



3 Continuum approaches 



The Id case Now let us focus on the continuum approaches describing 
the systems in terms of coarse grained velocity and density fields. As the 
simplest example, let us first investigate the Id SPP system described in the 
previous section. In [15] and [23] the following equations were obtained by the 
integration of the master equation of the microscopic dynamics: 

d tP = -v d x ( P U) + Dd 2 x p (12) 



and 



dt u = f(u) + p 2 d 2 x u + a MKM + c> 



(13) 



with p and U being the coarse-grained density and dimensionless velocity 
fileds, respectively. The coefficients vo, D, p and a are determined by the 
parameters of the microscopical rules of interaction. The function f(U) is 
antisymmetric with f(U) > for < U < 1 and f(U) < for U > 1. 
The noise ( has a zero mean and its standard deviation is proportional to 
p -1 / 2 . The coupling term in (13) with coefficient a is derived from the Tailor 
expansion of the local average velocity as 

lml , E-t U(x')p(x')dx' A : 



d 2 x U + 2 



(d x U)(d x p) 



+ ... (14) 



The nonlinear coupling term {d x U){d x p) / p is specific for such SPP systems: 
it is responsible for the slowing down (and eventually the "turning back") 
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of the particles under the influence of a larger number of particles moving 
oppositely. When two groups of particles move in the opposite direction, the 
density locally increases and the velocity decreases at the point they meet. 
Let us consider a particular case, when particles move from left to right and 
the velocity is locally decreasing while the density is increasing as x increases 
(particles are moving towards a "wall" formed between two oppositely moving 
groups). The term (d x u) is less, the term (d x p) is larger than zero in this case. 
Together they have a negative sign resulting in the slowing down of the local 
velocity. This is a consequence of the fact that there are more slower particles 
(in a given neighborhood) in the forward direction than faster particles coming 
from behind, so the average action experienced by a particle in the point x 
slows it down. 

For a = the dynamics of the velocity field U is independent of p, and Eq. (13) 
is equivalent to the time dependent Ginsburg-Landau <3> 4 model of spin chains, 
where domains of opposite magnetization develop at finite temperatures. As 
numerical simulations demonstrate [15,23], for large enough a and low noise 
the initial domain structure breaks and the system is organized into a single 
major group traveling in a spontaneously selected direction. This elimination 
of the domain structure, i.e., the ordering is strongly connected to the linear 
instability of the domains [15,23] for large enough a. 



Hydrodynamics in two and higher dimensions. To get an insight into 
the possible analytical treatment of SPP systems in higher dimensions, the 
Navier-Stokes equations can be applied for a fluid of SPPs, as described in 
[24]. The two basic equations governing the dynamics of the noiseless "self- 
propelled fluid" are the continuity equation 

d tP + V{pv) = (15) 



and the equation of motion 



d t v + (v-V)v = F(v) iTH «v> — u) Vp + u\/ 2 v, (16) 

T\ T 2 p* 

where p* is the effective density of the particles, p is the pressure, v is the 
kinematic viscosity, F{y) is the intrinsic driving force of biological origin and 
Ti, r 2 are time scales associated with velocity relaxation resulting from inter- 
action with the environment and the surrounding SPPs, respectively. Let us 
now go through the terms of Eq. 16. 
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The self-propulsion can be taken into account as a constant magnitude force 
acting in the direction of the velocity of the particles as 




T\ \V\ 



(17) 



where v o is the speed determined by the balance of the propulsion and friction 
forces, i.e., vo would be the speed of a homogeneous fluid. 

Taking Taylor expansion of the expression ^ ((v) — v) describing local velocity 
alignment yields 



similarly to the Id case (14). In the following we will consider cases where 
the density fluctuations are forced to be small. Hence the velocity-density 
coupling term in (18) is negligible, which is a major difference compared to 
the previously studied Id system. Thus, if the relative density changes are 
small, the velocity alignment can be incorporated with an effective viscosity 
coefficient v* into the viscous term of (16), 

In the simplest cases the pressure can be composed of an effective "hydro- 
static" pressure, and an externally applied pressure as 



where g* is a parameter related to the compressibility of the fluid. 

Combining (15), (16), (17) and (19) one gets the following final form for the 
equations of the SPP flow: 




— * 

— V 




(18) 



P = 9*P + Pext 



(19) 



d tP + V(pv) = 0, 



(20) 



and 



d t v + {v • V)v 



Vq V 



—v + u*V 2 v-g*Vp + — ^Vpcxt- 



(21) 



Ti \V 
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Fig. 6. (a) A numerically generated vortex for dimensionless system size R = 4.4 
and gi = 750. The length of the arrows is proportional to the local velocity, while 
their thickness is proportional to the density, (b) The measured (circles) and the 
theoretical (solid line) velocity profile for the vortex shown in (a). (After [24].) 

It is useful to introduce the characteristic length A = y/v*Ti and characteristic 
time T = ^ and write (20) and (21) in a dimensionless form of 

^ + VV) = 0, (22) 
r ± + $ . vV) = ^ - if + V'V - g[V P + g' 2 V Pcxt . (23) 



with t' = t/T, x' = x/X, v' = v/v and the V operator derivates with respect 
to f" — r/X. We also introduced the notations g[ = Trig* /A 2 and g' 2 = g[/g*- 
In the following we will drop the prime for simplicity. 

For certain simple geometries it is possible to obtain analytical stationary 
solutions [24] of the above equations supposing incompressibility (p = const.) 
and see that the minimal size of a vortex is in the order of A, i.e., one in 
dimensionless units. Therefore, if the dimensionless system size, R ^> 1 then 
initially many Kosterlitz-Thouless-like vortices are likely to be present in the 
system. 

Numerical solutions of (20) and (21) were also calculated in [24] with finite 
compressibility and p ext = 0. The only remaining dimensionless quantity char- 
acterizing the system is T\/T which relates the various relaxation times and 
v . The system can be considered overdamped if T\ T holds. Fig. 6 shows 
the stationary state for the equations in the high compressibility (gi = 750) 
and overdamped limit. The length and direction of the arrows show the ve- 
locity, while the thickness is proportional to the local density of the fluid. In 
Fig. 6b the radial velocity distribution is presented for the vortex shown in 
Fig. 6a together with the velocity profile obtained analytically in [24] . Rather 
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Fig. 7. Transient multiple vortex state at R = 25.3, g x = 128. (After [24].) 

good agreement is seen; the differences are due to the fact that the numerical 
system is not perfectly circular. In simulations with R ^> 1 multi- vortex states 
develop (Fig. 7), which transform into a single vortex after a long enough time. 

Thus, we demonstrated that - in contrast to the Id case - the ordering in 
2d is not due to the density- velocity coupling term in the expansion of (v): 
the viscosity and the internal driving force terms are sufficient to destabilize 
the originally present vortices and organize the system into a globally ordered 
phase. The stability of that ordered phase against a finite amount of fluctua- 
tions has been shown by Tu and Toner in [13,14] and discussed in this volume, 
too. 
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